The goals for this phase of your final project are:
Articulate an interesting research question based on a dataset you’d like to learn more about.
Develop a spatial database that contains potentially relevant explanatory variables that you’d like to explore in the context of that research question.
Demonstrate an understanding of the various workflow elements involved in designing and constructing a spatial database for subsequent visualization and analysis.
My Question: Where is the best place to go in New Zealand during the zombie apocalypse to avoid infection and survive? Specifically, I will look at 1. Where the most zombie attacks are happening 2. Are they happening in more urban or rural locations? 3. Highlight rural locations as areas of survival (the less people around, the less chance of encountering an infected person and becoming a zomibe) 4. Number of pharmacies in each city>>region to inform survivors of the best places to go where they will have regular access to medical supplies
I will answer this question using the following data:
library(here)
## here() starts at /home/michaela/R/Mini-Project-1-michaelagustafson
library(tidyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(terra)
## terra version 1.4.22
##
## Attaching package: 'terra'
## The following object is masked from 'package:dplyr':
##
## src
## The following object is masked from 'package:tidyr':
##
## extract
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidyverse)
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x terra::extract() masks tidyr::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::simplify() masks terra::simplify()
## x terra::src() masks dplyr::src()
library(nngeo)
library(data.table)
##
## Attaching package: 'data.table'
## The following object is masked from 'package:purrr':
##
## transpose
## The following objects are masked from 'package:terra':
##
## copy, shift
## The following objects are masked from 'package:dplyr':
##
## between, first, last
library(tmap)
## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sf
library(viridis)
## Loading required package: viridisLite
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
##
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
##
## project
# library(ggmap)
Data source: kaggle.com
zom <- read.csv("zombies.csv")
head(zom); dim(zom)
## zombieid zombie age sex rurality household water food medication tools
## 1 1 Human 18 Female Rural 1 0 Food Medication No tools
## 2 2 Human 18 Male Rural 3 24 Food Medication tools
## 3 3 Human 18 Male Rural 4 16 Food Medication No tools
## 4 4 Human 19 Male Rural 1 0 Food Medication tools
## 5 5 Human 19 Male Urban 1 0 Food Medication No tools
## 6 6 Human 19 Female Urban 1 0 Food Medication tools
## firstaid sanitation clothing documents
## 1 First aid supplies Sanitation Clothing <NA>
## 2 First aid supplies Sanitation Clothing <NA>
## 3 First aid supplies Sanitation Clothing <NA>
## 4 No first aid supplies Sanitation Clothing <NA>
## 5 First aid supplies Sanitation <NA> <NA>
## 6 First aid supplies Sanitation Clothing <NA>
## [1] 200 14
zombie.shp <- read_sf(here("New_Zealand_Zombie_Attack_Data.shp"))
plot(zombie.shp) #remember to use st_geometry to avoid plotting every column in the dataframe
Data source: https://datafinder.stats.govt.nz/layer/98765-regional-council-2019-clipped-generalised/
The regions of New Zealand are my ‘regions’ I will use to organize higher order variables.
# Find "Regional Council 2018 Clipped (generalised)"
# select the GeoPackage option in the "Vectors/tables" dropdown
# at https://datafinder.stats.govt.nz/data/ (requires registration)
# Save the result as:
unzip("statsnzregional-council-2019-clipped-generalised-SHP.zip")
nz.full.shp <- st_read("regional-council-2019-clipped-generalised.shp")
## Reading layer `regional-council-2019-clipped-generalised' from data source
## `/home/michaela/R/Mini-Project-1-michaelagustafson/regional-council-2019-clipped-generalised.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 17 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1089970 ymin: 4747987 xmax: 2470042 ymax: 6223156
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
plot(st_geometry(nz.full.shp))
# print(object.size(nz_full), units = "Kb") # 14407.2 Kb
# nz <- st_simplify(nz_full, dTolerance = 1000)
names(nz.full.shp)
## [1] "REGC2019_V" "REGC2019_1" "LAND_AREA_" "AREA_SQ_KM" "Shape_Leng"
## [6] "geometry"
nz.full.shp$REGC2019_1
## [1] "Northland Region" "Auckland Region"
## [3] "Waikato Region" "Bay of Plenty Region"
## [5] "Gisborne Region" "Hawke's Bay Region"
## [7] "Taranaki Region" "Manawatu-Wanganui Region"
## [9] "Wellington Region" "West Coast Region"
## [11] "Canterbury Region" "Otago Region"
## [13] "Southland Region" "Tasman Region"
## [15] "Nelson Region" "Marlborough Region"
## [17] "Area Outside Region"
nz.reg.shp = filter(nz.full.shp, REGC2019_1 != "Area Outside Region") %>%
select(Region = REGC2019_1, AreaSQKM = AREA_SQ_KM, geometry)
st_is_valid(nz.reg.shp)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE
st_crs(nz.reg.shp)
## Coordinate Reference System:
## User input: NZGD2000 / New Zealand Transverse Mercator 2000
## wkt:
## PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
## BASEGEOGCRS["NZGD2000",
## DATUM["New Zealand Geodetic Datum 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4167]],
## CONVERSION["New Zealand Transverse Mercator 2000",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",173,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",1600000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["New Zealand - onshore"],
## BBOX[-47.33,166.37,-34.1,178.63]],
## ID["EPSG",2193]]
Data source: https://www.stats.govt.nz/topics/census
I will use census data to determine the population density/sqkm in areas of the regions (urban vs rural) and overall population density/sqkm of a region. This information will be best used to determine good regions for individuals to go to be in less populated areas and have a better chance of avoiding other people and potential zombie infection.
# read in census data
nz.census <- read.csv("Census_Usually_resident_population_count_and_change_by_region_2006_2013_and_2018.csv")
# filter to get most recent census population (2018)
names(nz.census)
## [1] "X...Census.year" "Region" "Region.Code" "Measure"
## [5] "Value" "Value.Unit" "Value.Label" "Null.Reason"
nz.census <- filter(nz.census, X...Census.year == "2018") %>%
select(Year = X...Census.year, Region, Population = Value)
nz.census$Region = paste(nz.census$Region, "Region", sep = " ")
nz.census[8, 2] = "Manawatu-Wanganui Region"
You’ll want to leave yourself a few more “breadcrumbs”, I think. Anytime you are directly modifying the data (as you do above), it’s good to leave a note to remind you why you did
Data source: https://datafinder.stats.govt.nz/layer/98743-urban-rural-2019-clipped-generalised/
unzip("statsnzurban-rural-2019-clipped-generalised-SHP.zip")
urbrur <- read_sf("urban-rural-2019-clipped-generalised.shp")
plot(st_geometry(urbrur))
urban <- urbrur %>%
select(Name = UR2019_V_1, Urban = IUR2019__1, AreaSQKM = AREA_SQ_KM, geometry)
unique(urban$Urban)
## [1] "Inland water" "Large urban area" "Major urban area"
## [4] "Medium urban area" "Oceanic" "Rural other"
## [7] "Rural settlement" "Small urban area"
# neet to change columsn to say just urban or rural
urban$Urban[grep("urban", urban$Urban)] <- "Urban"
urban$Urban[grep("Rural", urban$Urban)] <- "Rural"
urban$Urban[grep("water", urban$Urban)] <- "Other"
urban$Urban[grep("Oceanic", urban$Urban)] <- "Other"
Data source: https://fyi.org.nz/request/10183-total-number-of-licenced-pharmacies-in-nz-2019
Does this still count as tabular data even though I turned it into spatial data??
# read in pharmacy csv with location data
nz.pharm <- read.csv(here("nz.pharm.csv"))
nz.pharm <- na.omit(nz.pharm)
nz.pharm.sf <- st_as_sf(nz.pharm, coords = c("lon", "lat"), crs = 4326)
st_crs(nz.pharm.sf)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
nz.pharm.proj <- st_transform(nz.pharm.sf, crs = st_crs(nz.reg.shp))
st_crs(nz.reg.shp) == st_crs(nz.pharm.proj)
## [1] TRUE
plot(st_geometry(nz.reg.shp))
plot(st_geometry(nz.pharm.proj), add = TRUE)
nz.pharm.proj$type <- "Pharmacy"
#looks alright
Spatial Join
Questions to Answer: 1. What is the percentage of urban and rural areas within each region? 2. What is the number of attacks in each urban type within each region? 3. What is the population density of each region (people/sqkm)? 4. What is the number of attacks per population in each region? (attacks/people)
# make sure CRS are the same
st_crs(zombie.shp) == st_crs(urban)
## [1] FALSE
zombie.shp.proj <- st_transform(zombie.shp, crs = st_crs(urban))
st_crs(zombie.shp.proj) == st_crs(urban)
## [1] TRUE
st_crs(urban) == st_crs(nz.reg.shp)
## [1] TRUE
# ended up using a nearest neighbor approach because I kept getting points 800 something rows when I started with 668 so something was getting duplicated:
#zombie.urban.join <- st_join(urban, zombie.shp.proj, join = st_nn, k = 1, maxdist = 50)
#zombie.urban.region.join <- st_join(zombie.urban.join, nz.reg.shp, join = st_within)
#str(zombie.urban.region.join)
#zombie.urban.region.ng <- st_drop_geometry(zombie.urban.region.join)
#zombie.urban.region.df <- zombie.urban.region.ng %>%
# select(Urban, AreaSQKM.x, Attacks, Region)
#zombie.urban.region.df <- as.data.frame(zombie.urban.region.df)
#zombie.urban.region.df$Region <- as.factor(zombie.urban.region.df$Region)
#zombie.urban.region.df$Attacks <- replace_na(zombie.urban.region.df$Attacks, 0)
#zombie.urban.region.df$Attacks <- as.numeric(zombie.urban.region.df$Attacks)
#summary.df <- zombie.urban.region.df %>%
#group_by(Region, Urban) %>%
#summarise(., sum(Attacks), sqkm_dis_total = sum(AreaSQKM.x))
#summary.df <- rename(summary.df, Attacks = "sum(Attacks)")
#summary.df$urban_attacks <- paste(summary.df$Urban, summary.df$Attacks, sep = ";")
#summary.df.wide <- spread(summary.df, Urban, sqkm_dis_total, "0")
#summary.df.wide <- spread(summary.df.wide, urban_attacks, Attacks)
#str(summary.df.wide)
#summary.df.wide$Other <- as.numeric(summary.df.wide$Other)
#summary.df.wide$Rural <- as.numeric(summary.df.wide$Rural)
#summary.df.wide$Urban <- as.numeric(summary.df.wide$Urban)
#summary.df.wide <- summary.df.wide[!is.na(summary.df.wide$Region),]
#summary.df.wide[is.na(summary.df.wide)] <- 0
#summary.df.wide <- summary.df.wide %>%
# group_by(Region) %>%
#summarise_if(., is.numeric, sum)
#summary.df.wide <- rename(summary.df.wide, Other_sqkm = "Other")
#summary.df.wide <- rename(summary.df.wide, Rural_sqkm = "Rural")
#summary.df.wide <- rename(summary.df.wide, Urban_sqkm = "Urban")
#
#colnames(summary.df.wide)[grepl('Other;',colnames(summary.df.wide))] <- 'Other_attacks'
#colnames(summary.df.wide)[grepl('Rural;',colnames(summary.df.wide))] <- 'Rural_attacks'
#colnames(summary.df.wide)[grepl('Urban;',colnames(summary.df.wide))] <- 'Urban_attacks'
#names(summary.df.wide) <- make.unique(names(summary.df.wide))
#tidy.summary.df <- summary.df.wide %>%
#rowwise() %>%
#mutate(Rural_Attack_sum = sum(across(starts_with("Rural_attacks")), na.rm = T))
#tidy.summary.df <- tidy.summary.df %>%
#rowwise() %>%
# mutate(Other_Attack_sum = sum(across(starts_with("Other_attacks")), na.rm = T))
#tidy.summary.df <- tidy.summary.df %>%
#rowwise() %>%
#mutate(Urban_Attack_sum = sum(across(starts_with("Urban_attacks")), na.rm = T))
# taking out Other Attack sum becuase there were none (they were in an unidentified region)
#tidy.summary.df <- tidy.summary.df %>%
#select(Region, Other_sqkm, Rural_sqkm, Urban_sqkm, Rural_Attack_sum, Urban_Attack_sum)
### NEED TO FIX REGIONAL SQKM
# Find the number of pharmacies per region
pharm.per.region <- st_join(nz.reg.shp, nz.pharm.proj, join = st_contains)
pharm.per.region.ng <- st_drop_geometry(pharm.per.region)
pharm.per.region.ng <- pharm.per.region.ng %>%
select(Region, Pharmacy_Nm)
# Count the number of pharmacies in each region
pharm.per.region.count <- pharm.per.region.ng %>%
group_by(Region) %>%
tally()
pharm.per.region.count <- rename(pharm.per.region.count, pharms_per_region = n)
# combine pharmacies per region with summary.df
#tidy.summary.df <- left_join(tidy.summary.df, pharm.per.region.count)
# add total population to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.census)
# add total area (sqkm) to final df
#tidy.summary.df <- left_join(tidy.summary.df, nz.reg.shp)
# calculate attacks/rural area and attacks per urban area per region
#tidy.summary.df <- tidy.summary.df %>%
#mutate(attacks_per_rural = Rural_Attack_sum/Rural_sqkm,
# attacks_per_urban = Urban_Attack_sum/Urban_sqkm,
# pop_density = Population/AreaSQKM,
# perc_rural = Rural_sqkm/AreaSQKM,
# perc_urb = Urban_sqkm/AreaSQKM,
# total_attacks = Rural_Attack_sum + Urban_Attack_sum)
<<<<<<< HEAD # Add elevation
rastlist1 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/", pattern = '.tif', full.names=TRUE)
rastlist1
## [1] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//aklnd_25r.tif"
## [2] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//ecape_25r.tif"
## [3] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//hksbay_25r.tif"
## [4] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//nthcape_25r.tif"
## [5] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//taranaki_25r.tif"
## [6] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//waikato_25r.tif"
## [7] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//well_25r.tif"
## [8] "/opt/data/MP/Gustafson_Michaela_MP/nzdem-north//whngrei_25r.tif"
rn1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/aklnd_25r.tif")
rn2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/ecape_25r.tif")
rn3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/hksbay_25r.tif")
rn4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/nthcape_25r.tif")
rn5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/taranaki_25r.tif")
rn6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/waikato_25r.tif")
rn7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/well_25r.tif")
rn8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem-north/whngrei_25r.tif")
rastlist2 <- list.files(path = "/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/", pattern = '.tif', all.files=TRUE, full.names=FALSE)
rastlist2
## [1] "chch25r.tif" "dunedin_25r.tif" "greymth_25r.tif" "invcgll_25r.tif"
## [5] "kaik_25r.tif" "mtcook_25r.tif" "nelson_25r.tif" "teanau_25r.tif"
## [9] "waitaki_25r.tif"
rs1 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/chch25r.tif")
rs2 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/dunedin_25r.tif")
rs3 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/greymth_25r.tif")
rs4 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/invcgll_25r.tif")
rs5 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/kaik_25r.tif")
rs6 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/mtcook_25r.tif")
rs7 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/nelson_25r.tif")
rs8 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/teanau_25r.tif")
rs9 <- rast("/opt/data/MP/Gustafson_Michaela_MP/nzdem_south/waitaki_25r.tif")
all.rast <- merge(rn1, rn2, rn3, rn4, rn5, rn6, rn7, rn8, rs1, rs2, rs3, rs4, rs5, rs6, rs7, rs8, rs9)
zombie.elev.proj <- st_transform(zombie.shp, crs = crs(all.rast))
plot(all.rast)
plot(zombie.elev.proj, add = TRUE)
## Warning in plot.sf(zombie.elev.proj, add = TRUE): ignoring all but the first
## attribute
# Technically should be calculating the mean elevation or range or something from a buffer surrounding the point, but for the sake of time (because I already tried it and I'd have to mess around with changing extents and a bunch of other stuff) and to actually practice what I'm supposed to be practicing - just going to the pull the elevation from the exact point of attack:
zomb.elev.proj.vect <- vect(zombie.elev.proj)
elev.at.attack <- extract(all.rast, zomb.elev.proj.vect, na.rm = TRUE)
elev.at.attack$ID <- zombie.shp$OBJECTID
elev.at.attack <- rename(elev.at.attack, OBJECTID = ID)
# add urban or rural
zombie.attacks.df <- st_join(zombie.shp.proj, urban, join = st_within)
colnames(zombie.attacks.df)
## [1] "OBJECTID" "OBJECTID_1" "Day" "Attacks" "Time"
## [6] "Longitude" "Latitude" "geometry" "Name" "Urban"
## [11] "AreaSQKM"
zombie.attacks.df <- zombie.attacks.df[,c(1, 4, 6:8, 10)]
# add elevation
zombie.attacks.df <- left_join(zombie.attacks.df, elev.at.attack)
## Joining, by = "OBJECTID"
zombie.attacks.df <- rename(zombie.attacks.df, elev = aklnd_25r)
# add region
zombie.attacks.df <- st_join(zombie.attacks.df, nz.reg.shp, join = st_within)
# add number of pharmacies per region
zombie.attacks.df <- left_join(zombie.attacks.df, pharm.per.region.count)
## Joining, by = "Region"
# add census data
zombie.attacks.df <- left_join(zombie.attacks.df, nz.census)
## Joining, by = "Region"
# mutate and calculate population density per region
zombie.attacks.df <- zombie.attacks.df %>%
mutate(pop_density = Population/AreaSQKM)
#####
nz.reg.shp.join <- left_join(nz.reg.shp, nz.census)
## Joining, by = "Region"
nz.reg.shp.join <- nz.reg.shp.join %>%
mutate(pop_density = Population/AreaSQKM)
======= you didn’t have any raster extractions here and the analysis is not reproducible because a) you didn’t save the nz_pharm object, b) you didn’t provide all of the necessary info to recreate it. I would have liked to see a little more effort and annotation. If I knew more about what you were trying to do, I might actually be able to fix the parts of this that wont run.
save.image(here("workspaces/MiniProj_Workspace.RData"))
library(tidyverse)
library(pander)
##
## Attaching package: 'pander'
## The following object is masked from 'package:terra':
##
## wrap
library(sf)
library(units)
## udunits database from /usr/share/xml/udunits/udunits2.xml
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
##
## Attaching package: 'ggmap'
## The following object is masked from 'package:terra':
##
## inset
library(cartogram)
##
## Attaching package: 'cartogram'
## The following object is masked from 'package:terra':
##
## cartogram
library(patchwork)
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:terra':
##
## area
library(tmap)
library(viridis)
library(tigris)
## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
library(ggspatial)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggmap':
##
## wind
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
zombie.attacks.df <- st_as_sf(zombie.attacks.df)
?tm_polygons
tm_shape(nz.reg.shp.join) +
tm_polygons(col = "Region", palette = "viridis") +
tm_legend(outside = TRUE) +
tm_bubbles(size = "Population", col = "pop_density", border.col = "black")
## Some legend labels were too wide. These labels have been resized to 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## The legend is too narrow to place all symbol sizes.
# required libraries
library(leaflet, quietly = T, warn.conflicts = F)
# start basemap
map <- leaflet() %>%
# add ocean basemap
addProviderTiles(providers$Esri.OceanBasemap) %>%
# add another layer with place names
addProviderTiles(providers$Hydda.RoadsAndLabels, group = 'Place names') %>%
# add graticules from a NOAA webserver
addWMSTiles(
"https://gis.ngdc.noaa.gov/arcgis/services/graticule/MapServer/WMSServer/",
layers = c("1-degree grid", "5-degree grid"),
options = WMSTileOptions(format = "image/png8", transparent = TRUE),
attribution = NULL,group = 'Graticules') %>%
# focus map in a certain area / zoom level
setView(lng = 174, lat = -41, zoom = 5)%>%
# add layers control
addLayersControl(overlayGroups = c('Place names',
'Graticules',
'Points',
'Lines',
'Polygons'),
options = layersControlOptions(collapsed = FALSE),
position = 'topright') %>%
# list groups to hide on startup
hideGroup(c('Place names'))
# show map
map
### Add data
nz.shp.reproj <- st_transform(nz.reg.shp.join, crs = st_crs(zombie.shp))
library(htmltools)
##
## Attaching package: 'htmltools'
## The following object is masked from 'package:pander':
##
## p
labs <- as.list(zombie.shp$Attacks)
labs2 <- as.list(nz.shp.reproj$pop_density)
## Color setup
colpal <- colorNumeric(palette = "magma", domain=nz.shp.reproj[['pop_density']], n=10)
colorData <- nz.shp.reproj[["pop_density"]]
map <- map %>%
addTiles() %>%
addPolygons(data = nz.shp.reproj, label = lapply(labs2, HTML), color = ~colpal(colorData))
map <- map %>%
addCircleMarkers(data = zombie.shp,
weight = 0.5,
label = lapply(labs, HTML),
col = 'black',
fillColor = 'darkslategrey',
radius = 4,
fillOpacity = 0.9,
stroke = T)
map